Manuscript title: Remdesivir prophylaxis but not later treatment limits measles-induced “immune amnesia” and measles antibody responses in macaques
Abstract: Despite the availability of a safe and effective vaccine, there is still no licensed antiviral therapy available for measles. We therefore evaluated the efficacy of a broad-spectrum antiviral, remdesevir, in wild-type measles virus (WT MeV)-infected RMs, and examines plasma antibody reactivity to viral epitopes via VirScan, a phage-display immunoprecipitation and sequencing (PhIP-Seq) technology. Here, VirScan analysis was performed on the longitudinal plasma samples from four RM treatment groups with both the RM and human virome library.
library(dplyr)
library(stringr)
library(tidyr)
library(ggplot2)
library(readr)
library(RColorBrewer)
library(gridExtra)
library(ggh4x)
library(ggpubr)
library(gt)
VirScan raw datasets were first cleaned to remove control beads information, and transformed to four working datasets:
Column explanation:
Treatment: Three groups of RMs were infected with WT MeV at D0, with the first group being given PBS as treatment (WTMeV), second group being given daily remdesevir infusion as post-exposure prohylasxis from D3 to D11 (RDV_PEP) and the third group being given daily remdesevir infusion late in infection from D11 to D22 (RDV_Late). The forth group was included as a sham control group (Sham_control).
DPI: Days post infusion (D0, D14_15, D28, D84_89, D112_168_176)
Macaques_ID: Identification number of each RM
# Load cleaned dataset
# VARscore_1 <- readRDS("VARscore_1.RDS")
# Merged_hfc_all <- readRDS("Merged_hfc_all.RDS")
# Total_peptide_hits_all <- readRDS("Total_peptide_hits_all.RDS")
# Total_number_of_peptide_tiles_per_pathogen_all_table
Total_number_of_peptide_tiles_per_pathogen_all <- Merged_hfc_all %>%
group_by(taxon_species, pep_aa, pos_start, pos_end, pep_id)%>%
summarize()
## `summarise()` has grouped output by 'taxon_species', 'pep_aa', 'pos_start',
## 'pos_end'. You can override using the `.groups` argument.
Total_number_of_peptide_tiles_per_pathogen_all_table <- as.data.frame(table(Total_number_of_peptide_tiles_per_pathogen_all$taxon_species))%>%
rename(taxon_species= Var1, Number_of_Peptide_Tiles = Freq)%>%
filter(!row_number() %in% c(1))
# Percentage peptide hits (new metric)
Percentage_peptide_hits_all <- Total_peptide_hits_all %>%
inner_join(Total_number_of_peptide_tiles_per_pathogen_all_table, by = "taxon_species") %>%
filter(Number_of_Peptide_Tiles > 39) %>%
mutate(pph = total_peptide_hits/Number_of_Peptide_Tiles*100)
RDV_groups <- c("WTMeV","RDV_PEP", "RDV_Late", "Sham_control")
color_pattern_for_RDV_groups <- c("red","pink", "#6600CC", "#006699")
VARscore_per_group_heatmap_function <- function (virus_group){
VARscore_subset_mean <- VARscore_1 %>%
filter(taxon_species %in% virus_group)%>%
group_by(Treatment, taxon_species, DPI)%>%
summarise(average_VAR = mean(VAR))
print(ggplot(VARscore_subset_mean, aes(x = DPI, y = taxon_species))+
geom_raster(aes(fill = average_VAR))+
scale_fill_gradient2(low="white", mid="white", high="darkblue", midpoint=0) +
labs(x = "DPI", y= "viral species", title = "Average VARscore for each RM group") +
facet_wrap(~Treatment, ncol = 4, scales = "fixed") +
theme_classic(base_size=14)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))
}
VARscore_per_RM_heatmap_function <- function (datatable, virus_group, treatment_group, color_pattern){
VARscore_per_RM_heatmap_graph_function <- function (datatable, virus_group, treatment_group, color_pattern) {
VARscore_subset <- VARscore_1 %>%
filter(taxon_species %in% virus_group) %>%
filter(Treatment == treatment_group)
print(ggplot(VARscore_subset, aes(x=Macaques_ID , y = taxon_species))+
geom_raster(aes(fill = VAR))+
scale_fill_gradient2(low ="white", high= color_pattern) +
labs(x = "Monkey ID", y= "viral species", title = paste0("Average VARscore in the ", treatment_group, " group")) +
facet_wrap( ~DPI, ncol = 5, nrow = 3, scales = "fixed")+
theme_classic(base_size=14)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))
}
VARscore_per_RM_heatmap_graph_function(datatable, virus_group, treatment_group[1], color_pattern[1])
VARscore_per_RM_heatmap_graph_function(datatable, virus_group, treatment_group[2],color_pattern[2])
VARscore_per_RM_heatmap_graph_function(datatable, virus_group, treatment_group[3], color_pattern[3])
VARscore_per_RM_heatmap_graph_function(datatable, virus_group, treatment_group[4], color_pattern[4])
}
VARscore_per_group_heatmap_function(top_viral_species_human_lib_2)
VARscore_per_group_heatmap_function(NHP_virus)
VARscore_per_RM_heatmap_function(VARscore_1, top_viral_species_human_lib_2, RDV_groups, color_pattern_for_RDV_groups)
VARscore_per_RM_heatmap_function(VARscore_1, NHP_virus, RDV_groups, color_pattern_for_RDV_groups)
VAR_graph_function_1 <- function (datatable) {
ggplot(datatable, aes(x = DPI, y = VAR, group = Macaques_ID))+
geom_line(aes(color = taxon_species), linewidth = 1, colour="lightgrey")+
labs(x = "DPI", y= "VAR Score") +
facet_wrap(~Treatment, ncol = 4, scales = "fixed") +
theme_classic(base_size=10)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
stat_summary(aes(group=Treatment), fun=mean, geom="line", colour="red", linewidth = 2)
}
VAR_graph_function_2 <- function (datatable_0) {
ggplot(datatable_0, aes(x=DPI, y=VAR, color = Treatment)) +
geom_point() +
theme_classic(base_size=10)+
geom_line(aes(group = Treatment), stat = "summary" , fun = mean) +
labs(x= "DPI" , y="Antibody reactivity (VAR scores)")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0.2))
}
VAR_graph_function_3 <- function (datatable_00) {
ggplot(datatable_00, aes(x=DPI, y=VAR, color = Treatment)) +
geom_boxplot()+
geom_point(position=position_jitterdodge(jitter.width = 0, jitter.height = 0))+
facet_wrap(~Treatment, ncol =4)+
labs(x= "DPI" , y="Antibody reactivity (VAR scores)")+
theme_classic(base_size=10)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0.2))+
theme(legend.position = "none")+
stat_compare_means(method = "kruskal.test", label.x = 1.5, label.y = -1, color="black", size = 2)
}
VAR_overall_function <- function (i){
plot_list_1 <<- list()
plot_list_2 <<- list()
plot_list_3 <<- list()
for (k in i){
datatable_1 <- VARscore_1 %>%
filter(taxon_species == k)
plot1 <- VAR_graph_function_1(datatable_1)+ labs(title = paste0("Changes in VAR scores of \n", k))
plot2 <- VAR_graph_function_2(datatable_1)+ labs(title = paste0("Antibody reactivity of \n", k))
plot3 <- VAR_graph_function_3(datatable_1)+ labs(title = paste0("Antibody reactivity of \n", k))
plot_list_1 [[k]] <<- plot1
plot_list_2 [[k]] <<- plot2
plot_list_3 [[k]] <<- plot3
}
}
VAR_overall_function (top_viral_species_human_lib_2)
grid.arrange(grobs = plot_list_1, ncol=3, nrow = 6)
grid.arrange(grobs = plot_list_2, ncol=3, nrow = 6)
grid.arrange(grobs = plot_list_3, ncol=3, nrow = 6)
VAR_overall_function(NHP_virus)
grid.arrange(grobs = plot_list_1, ncol=3, nrow = 2)
grid.arrange(grobs = plot_list_2, ncol=3, nrow = 2)
grid.arrange(grobs = plot_list_3, ncol=3, nrow = 2)
# MeV VARscore head on comparision per timepoint
ggplot(subset(VARscore_1, Treatment != "Sham_control" & taxon_species == "Measles virus"), aes(x = factor (Treatment, levels = c("WTMeV", "RDV_PEP", "RDV_Late")), y = VAR, fill = Treatment)) +
geom_boxplot()+
scale_fill_manual(values = c("#F8766D" , "#7CAE00" ,"#00BFC4")) +
stat_compare_means(method = "wilcox.test", comparison = list (c("WTMeV", "RDV_PEP"), c("WTMeV", "RDV_Late"), c("RDV_PEP", "RDV_Late")), size = 4, label = "p.signif") +
theme_minimal(base_size = 16) +
facet_wrap(~DPI, ncol =5)+
labs(title = "MeV VAR scores per timepoint per group",
x = "RM Group",
y = "VARScores",
fill = "RM Group")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0.2))
# MeV peptide hits head on comparision per timepoint
ggplot(subset(Total_peptide_hits_all, Treatment != "Sham_control" & taxon_species == "Measles virus"), aes(x = factor (Treatment, levels = c("WTMeV", "RDV_PEP", "RDV_Late")), y = total_peptide_hits, fill = Treatment)) +
geom_boxplot()+
scale_fill_manual(values = c("#F8766D" , "#7CAE00" ,"#00BFC4")) +
stat_compare_means(method = "wilcox.test", comparison = list (c("WTMeV", "RDV_PEP"), c("WTMeV", "RDV_Late"), c("RDV_PEP", "RDV_Late")), size = 4, label = "p.signif") +
theme_minimal(base_size = 16) +
facet_wrap(~DPI, ncol =5)+
labs(title = "MeV peptide hits per timepoint per group",
x = "RM Group",
y = "Total peptide hits",
fill = "RM Group")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0.2))
RM_id <- unique(Total_peptide_hits_all$Macaques_ID)
WTMeV_RM_id<- c("104F", "108F", "42F" , "43F", "80F", "92F")
RDV_PEP_RM_id <- c("35F", "41F", "64E", "79F", "83F")
RDV_Late_RM_id <- c("57F","88F", "90F")
Sham_control_RM_id <- c("117H", "177H", "83H", "84H")
color_pattern_for_RDV_groups <- c("#F8766D" , "#7CAE00" ,"#00BFC4","#C77CFF")
df_data <- data.frame()
VAR_numeric_change_over_time_per_RM <- function (datatable, color_pattern){
ggplot(datatable, aes(x =DPI, y = numeric_change, fill = Treatment))+
geom_boxplot(fill = color_pattern)+
geom_point(color = "grey")+
geom_line(aes(group = taxon_species), stat = "summary" , fun = mean, color = "grey") +
theme_classic(base_size=20)+
ylim(-2,1.5)+
theme(axis.text=element_text(size=20),
axis.title =element_text(size=20))+
labs(x= "DPI" , y="numeric change")+
theme(legend.position="none")
}
VAR_numeric_change_over_time_per_RM_overall_function <- function (RM_id_all, color_pattern, timepoint, timepoint_2){
plot_list_1 <<- list()
for (k in RM_id_all){
top_viral_species_based_on_VAR_variable <- VARscore_1 %>%
filter (Macaques_ID == k & VAR >1 & DPI == "d0")
top_viral_species_based_on_VAR_variable <- unique(top_viral_species_based_on_VAR_variable$taxon_species)
VAR_numeric_change_over_time_per_RM_dataset <- VARscore_1 %>%
filter(Macaques_ID == k & taxon_species %in% top_viral_species_based_on_VAR_variable & DPI %in% c("d0", timepoint)) %>%
filter(taxon_species != "Measles virus")%>%
select(-total_peps)%>%
ungroup() %>%
pivot_wider(names_from = DPI, values_from = VAR)%>%
mutate(Numeric_VAR_score_changes_from_D0 = {{timepoint_2}} - d0)%>%
mutate(d0 = 0)%>%
pivot_longer(col = c("d0", "Numeric_VAR_score_changes_from_D0"), names_to= "DPI", values_to= "numeric_change")
plot1 <- VAR_numeric_change_over_time_per_RM(VAR_numeric_change_over_time_per_RM_dataset, color_pattern)+
labs (title = paste0( k, "\nAvg numeric change = ", round(mean(VAR_numeric_change_over_time_per_RM_dataset$numeric_change)*2, 1),"\n", "Total viral species = ", length(unique(VAR_numeric_change_over_time_per_RM_dataset$taxon_species))))+
scale_x_discrete(labels = c("d0", timepoint))
plot_list_1 [[k]] <<- plot1
df_data <<- rbind(df_data, round(mean(VAR_numeric_change_over_time_per_RM_dataset$numeric_change)*2, 1))
}
}
df_data <- data.frame()
VAR_numeric_change_over_time_per_RM_overall_function (WTMeV_RM_id, "#F8766D","d14_15", d14_15 )
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (RDV_PEP_RM_id, "#7CAE00","d14_15", d14_15)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (RDV_Late_RM_id, "#00BFC4","d14_15", d14_15)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (Sham_control_RM_id, "#C77CFF", "d14_15", d14_15)
grid.arrange(grobs = plot_list_1, ncol = 6)
df_data_1 <- df_data
df_data <- data.frame()
VAR_numeric_change_over_time_per_RM_overall_function (WTMeV_RM_id, "#F8766D","d28", d28 )
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (RDV_PEP_RM_id, "#7CAE00","d28", d28)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (RDV_Late_RM_id, "#00BFC4","d28", d28)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (Sham_control_RM_id, "#C77CFF", "d28", d28)
grid.arrange(grobs = plot_list_1, ncol = 6)
df_data_2 <- df_data
df_data <- data.frame()
VAR_numeric_change_over_time_per_RM_overall_function (WTMeV_RM_id[1:5], "#F8766D","d84_89", d84_89 )
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (RDV_PEP_RM_id, "#7CAE00","d84_89", d84_89)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (RDV_Late_RM_id, "#00BFC4","d84_89", d84_89)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (Sham_control_RM_id, "#C77CFF", "d84_89", d84_89)
grid.arrange(grobs = plot_list_1, ncol = 6)
df_data_3 <- df_data
df_data <- data.frame()
VAR_numeric_change_over_time_per_RM_overall_function (WTMeV_RM_id, "#F8766D", "d112_168_176", d112_168_176)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (RDV_PEP_RM_id, "#7CAE00", "d112_168_176", d112_168_176)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (RDV_Late_RM_id, "#00BFC4", "d112_168_176", d112_168_176)
grid.arrange(grobs = plot_list_1, ncol = 6)
VAR_numeric_change_over_time_per_RM_overall_function (Sham_control_RM_id, "#C77CFF", "d112_168_176", d112_168_176)
grid.arrange(grobs = plot_list_1, ncol = 6)
df_data_4 <- df_data
df_data_1 <- df_data_1 %>%
mutate(group = "d0 to d14_15")%>%
mutate(RM_ID = c("104F", "108F", "42F" , "43F", "80F", "92F", "35F", "41F", "64E", "79F", "83F", "57F","88F", "90F", "117H", "177H", "83H", "84H")) %>%
mutate(RM_group = c("WTMeV", "WTMeV", "WTMeV", "WTMeV", "WTMeV", "WTMeV", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_Late", "RDV_Late", "RDV_Late", "Sham_control", "Sham_control", "Sham_control", "Sham_control"))%>%
rename (average_VAR_numeric_change = 1)%>%
filter(RM_ID !="177H" )
df_data_2 <- df_data_2 %>%
mutate(group = "d0 to d28")%>%
mutate(RM_ID = c("104F", "108F", "42F" , "43F", "80F", "92F", "35F", "41F", "64E", "79F", "83F", "57F","88F", "90F", "117H", "177H", "83H", "84H")) %>%
mutate(RM_group = c("WTMeV", "WTMeV", "WTMeV", "WTMeV", "WTMeV", "WTMeV", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_Late", "RDV_Late", "RDV_Late", "Sham_control", "Sham_control", "Sham_control", "Sham_control"))%>%
rename (average_VAR_numeric_change = 1)%>%
filter(RM_ID !="177H" )
df_data_3 <- df_data_3 %>%
mutate(group = "d0 to d84_89")%>%
mutate(RM_ID = c("104F", "108F", "42F" , "43F", "80F", "35F", "41F", "64E", "79F", "83F", "57F","88F", "90F", "117H", "177H", "83H", "84H")) %>%
mutate(RM_group = c("WTMeV", "WTMeV", "WTMeV", "WTMeV", "WTMeV", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_Late", "RDV_Late", "RDV_Late", "Sham_control", "Sham_control", "Sham_control", "Sham_control"))%>%
rename (average_VAR_numeric_change = 1)%>%
filter(RM_ID !="177H" )
df_data_4 <- df_data_4 %>%
mutate(group = "d0 to d112_168_176")%>%
mutate(RM_ID = c("104F", "108F", "42F" , "43F", "80F", "92F", "35F", "41F", "64E", "79F", "83F", "57F","88F", "90F", "117H", "177H", "83H", "84H")) %>%
mutate(RM_group = c("WTMeV", "WTMeV", "WTMeV", "WTMeV", "WTMeV", "WTMeV", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_PEP", "RDV_Late", "RDV_Late", "RDV_Late", "Sham_control", "Sham_control", "Sham_control", "Sham_control"))%>%
rename (average_VAR_numeric_change = 1)%>%
filter(RM_ID !="177H" )
WT_PEP_Late_group_comparison <- function (dataset) {
ggplot(dataset, aes(x = factor (RM_group, levels = c("WTMeV", "RDV_PEP", "RDV_Late", "Sham_control")), y = as.numeric(average_VAR_numeric_change), fill = RM_group)) +
geom_boxplot(fill = c("#F8766D" , "#7CAE00" ,"#00BFC4","#C77CFF")) +
stat_compare_means(method = "wilcox.test", comparison = list (c("WTMeV", "RDV_PEP"), c("WTMeV", "RDV_Late"), c("RDV_PEP", "RDV_Late"), c("WTMeV", "Sham_control"), c("RDV_PEP", "Sham_control"), c("RDV_Late", "Sham_control")), size = 4, label = "p.signif") +
theme_minimal(base_size = 14) +
labs(title = paste0(unique(dataset$group)),
x = "RM Group",
y = "Average Numeric Change",
fill = "RM Group") +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
}
WT_PEP_Late_group_comparison(df_data_1)
WT_PEP_Late_group_comparison(df_data_2)
WT_PEP_Late_group_comparison(df_data_3)
WT_PEP_Late_group_comparison(df_data_4)
df_data_all <- rbind (df_data_1, df_data_2, df_data_3, df_data_4)
ggplot(df_data_all, aes(x = factor (group , levels = c("d0 to d14_15", "d0 to d28", "d0 to d84_89", "d0 to d112_168_176")), y = as.numeric(average_VAR_numeric_change), fill = factor (group , levels = c("d0 to d14_15", "d0 to d28", "d0 to d84_89", "d0 to d112_168_176")))) +
geom_boxplot() +
stat_compare_means(method = "kruskal.test", size = 3, paired = TRUE) +
facet_wrap (~factor (RM_group, levels = c("WTMeV", "RDV_PEP", "RDV_Late", "Sham_control")))+
theme_minimal(base_size = 16) +
labs(x = "RM Group",
y = "Average Numeric Change",
fill = "RM Group")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
ggplot(df_data_all, aes(x = factor (group , levels = c("d0 to d14_15", "d0 to d28", "d0 to d84_89", "d0 to d112_168_176")), y = as.numeric(average_VAR_numeric_change), color = factor (RM_group, levels = c("WTMeV", "RDV_PEP", "RDV_Late", "Sham_control")))) +
geom_boxplot() +
scale_color_manual(values = c("#F8766D" , "#7CAE00" ,"#00BFC4","#C77CFF"))+
stat_compare_means(method = "kruskal.test", size = 3, paired = TRUE) +
facet_wrap (~factor (RM_group, levels = c("WTMeV", "RDV_PEP", "RDV_Late", "Sham_control")), ncol = 4)+
theme_minimal(base_size = 14) +
labs(x = "DPI",
y = "Average Numeric Change",
color = "RM Group")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
PPH_per_group_heatmap_function <- function (datatable, virus_group, color_pattern) {
PPH_variable <-datatable%>%
filter(taxon_species %in% virus_group) %>%
group_by(taxon_species,DPI,Treatment)%>%
mutate(average_pph = mean(pph))
print(ggplot(PPH_variable, aes(x=DPI , y = taxon_species))+
geom_raster(aes(fill = average_pph))+
scale_fill_gradient2(low ="white", high= color_pattern) +
labs(x = "DPI", y= "viral species", title = "Average pph for each RM group") +
facet_wrap( ~Treatment, ncol = 5, nrow = 1, scales = "fixed")+
theme_classic(base_size=14)+
theme(legend.position="right")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))
}
PPH_per_RM_heatmap_function <- function (datatable, virus_group, treatment_group, color_pattern){
PPH_per_RM_graph_function <- function (datatable, virus_group, treatment_group, color_pattern) {
PPH_variable <-datatable%>%
filter(taxon_species %in% virus_group) %>%
filter(Treatment == treatment_group)
print(ggplot(PPH_variable, aes(x=Macaques_ID , y = taxon_species))+
geom_raster(aes(fill = pph))+
scale_fill_gradient2(low ="white", high= color_pattern) +
labs(x = "Monkey ID", y= "viral species", title = paste0("Percentage peptide hits in the ", treatment_group, " group")) +
facet_wrap( ~DPI, ncol = 5, nrow = 3, scales = "fixed")+
theme_classic(base_size=14)+
theme(legend.position="right")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))
}
PPH_per_RM_graph_function(datatable, virus_group, treatment_group[1], color_pattern[1])
PPH_per_RM_graph_function(datatable, virus_group, treatment_group[2],color_pattern[2])
PPH_per_RM_graph_function(datatable, virus_group, treatment_group[3], color_pattern[3])
PPH_per_RM_graph_function(datatable, virus_group, treatment_group[4], color_pattern[4])
}
RDV_groups <- c("WTMeV","RDV_PEP", "RDV_Late", "Sham_control")
color_pattern_for_RDV_groups <- c("red","pink", "#6600CC", "#006699")
PPH_per_group_heatmap_function(Percentage_peptide_hits_all, top_viral_species_human_lib_10, "springgreen4")
PPH_per_group_heatmap_function(Percentage_peptide_hits_all, NHP_virus, "springgreen4")
PPH_per_RM_heatmap_function(Percentage_peptide_hits_all, top_viral_species_human_lib_10, RDV_groups, color_pattern_for_RDV_groups)
PPH_per_RM_heatmap_function(Percentage_peptide_hits_all, NHP_virus, RDV_groups, color_pattern_for_RDV_groups)
PPH_graph_function <- function (datatable) {
ggplot(datatable, aes(x=DPI, y=pph, color = Treatment)) +
geom_point() +
geom_line(aes(group = Treatment), stat = "summary" , fun = mean) +
theme_classic(base_size=12)+
labs(x= "DPI" , y="Percentage peptide hits")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0.2))
}
PPH_graph_function_2 <- function (datatable) {
ggplot(datatable, aes(x=DPI, y=pph, color = Treatment)) +
geom_boxplot()+
geom_point(position=position_jitterdodge(jitter.width = 0, jitter.height = 0))+
facet_wrap(~Treatment, ncol=4)+
labs(x= "DPI" , y="Percentage peptide hits")+
theme_classic(base_size=12)+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=0.2))+
stat_compare_means(method = "kruskal.test", label.x = 1.5, label.y = -2, color="black", size = 2.5)
}
PPH_overall_function <- function (virus_group){
plot_list_1 <<- list()
plot_list_2 <<- list()
for (k in virus_group){
datatable<- Percentage_peptide_hits_all %>%
filter(taxon_species == k)
plot1 <- PPH_graph_function(datatable)+ labs(title = paste0("PPH of \n", k))
plot2 <- PPH_graph_function_2(datatable)+ labs(title = paste0("PPH of \n", k))
plot_list_1 [[k]] <<- plot1
plot_list_2 [[k]] <<- plot2
}
}
PPH_overall_function (top_viral_species_human_lib_10)
grid.arrange(grobs = plot_list_1, ncol = 3, nrow = 8)
grid.arrange(grobs = plot_list_2, ncol = 2, nrow = 12)
PPH_overall_function (NHP_virus)
grid.arrange(grobs = plot_list_1, ncol=3, nrow = 2)
grid.arrange(grobs = plot_list_2, ncol=2, nrow = 3)
RDV_groups <- c("WTMeV","RDV_PEP", "RDV_Late", "Sham_control")
color_pattern_for_RDV_groups <- c("red","pink", "#6600CC", "#006699")
Viral_peptide_hits_datatable_function <- function(datatable, virus) {
peptide_hits <-datatable %>%
filter(taxon_species == virus)%>%
filter_at(vars(gene_symbol), all_vars(!is.na(.)))%>%
mutate(hit_or_not = ifelse(hfc > 1, 1,0))%>%
group_by(Macaques_ID, gene_symbol,DPI,Treatment)%>%
mutate(total_peptide_hits= sum(hit_or_not)) %>%
group_by(gene_symbol,DPI,Treatment)%>%
mutate(average_peptide_hits = mean(total_peptide_hits))
return (peptide_hits)
}
Average_peptide_hits_per_protein_function <- function (datatable, virus, color_pattern){
print(ggplot(datatable, aes(x=DPI , y = gene_symbol))+
geom_raster(aes(fill = average_peptide_hits))+
scale_fill_gradient2(low ="white", high=color_pattern,
midpoint=0) +
labs(x = "DPI", y= "Protein", title = paste0(virus," epitope mapping"), fill = "Average peptide hit per protein") +
facet_wrap( ~Treatment, ncol = 4, scales = "fixed")+
theme_classic(base_size=15)+
theme(legend.position="bottom")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))
}
Peptide_hits_per_protein_function <- function (datatable, virus, treatment_group, color_pattern) {
Peptide_hits_variable <-datatable %>%
filter(Treatment==treatment_group)
print(ggplot(Peptide_hits_variable, aes(x=Macaques_ID , y = gene_symbol))+
geom_raster(aes(fill = total_peptide_hits))+
scale_fill_gradient2(low ="white", high=color_pattern,
midpoint=0) + #change to 10 or 12
labs(x = "RM ID", y= paste0 (virus, " proteins"), title = paste0("Total peptide hits in ",treatment_group), fill = "Total peptide hit \nper protein") +
facet_wrap(Treatment ~ DPI, ncol = 5, nrow = 1, scales = "fixed")+
theme_classic(base_size=15)+
theme(legend.position="right")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))
}
Peptide_hits_per_protein_function_all <- function (datatable, virus, treatment_group, color_pattern) {
Peptide_hits_variable <-datatable %>%
filter(Treatment==treatment_group)
print(ggplot(Peptide_hits_variable, aes(x=Macaques_ID , y = gene_symbol))+
geom_raster(aes(fill = total_peptide_hits))+
scale_fill_gradient2(low ="white", high=color_pattern,
midpoint=0) + #change to 10 or 12
labs(x = "RM ID", y= paste0 (virus, " proteins"), title = paste0("Total peptide hits in ",treatment_group), fill = "Total peptide hit \nper protein") +
facet_wrap(Treatment ~ DPI, ncol = 5, nrow = 2, scales = "fixed")+
theme_classic(base_size=12)+
theme(legend.position="right")+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))
}
Overall_peptide_hits_per_protein_function<- function (datatable, virus, treatment_group, color_pattern){
Peptide_hits_per_protein_function(datatable, virus, treatment_group[1], color_pattern[1])
Peptide_hits_per_protein_function(datatable, virus, treatment_group[2],color_pattern[2])
Peptide_hits_per_protein_function(datatable, virus, treatment_group[3], color_pattern[3])
Peptide_hits_per_protein_function(datatable, virus, treatment_group[4], color_pattern[4])
}
Overall_peptide_hits_per_protein_function_all<- function (datatable, virus, treatment_group, color_pattern){
Peptide_hits_per_protein_function_all(datatable, virus, treatment_group[1], color_pattern[1])
Peptide_hits_per_protein_function_all(datatable, virus, treatment_group[2],color_pattern[2])
Peptide_hits_per_protein_function_all(datatable, virus, treatment_group[3], color_pattern[3])
Peptide_hits_per_protein_function_all(datatable, virus, treatment_group[4], color_pattern[4])
}
Measles_peptide_hits <- Viral_peptide_hits_datatable_function(Merged_hfc_all, "Measles virus") %>% #type in any virus of interest to visualize peptide hit per protein
mutate (gene_symbol = if_else(gene_symbol == "P/V", "P/V/C", gene_symbol))%>% #change problematic gene_symbol naming
mutate (gene_symbol = if_else(gene_symbol == "P", "P/V/C", gene_symbol))
Average_peptide_hits_per_protein_function(Measles_peptide_hits, "Measles", "springgreen4")
Overall_peptide_hits_per_protein_function_all(Measles_peptide_hits,"Measles", RDV_groups, color_pattern_for_RDV_groups)
#RDV facet designs
design <- "
ABCDEF
GHIJK#
LMN###
OPQR##
"
Epitopes_list_and_dataset_function <- function (datatable, virus, hfc_value){
Epitopes_list <- datatable %>%
filter(taxon_species == virus) %>%
# mutate (gene_symbol = if_else(gene_symbol == "P/V", "P/V/C", gene_symbol))%>%
# mutate (gene_symbol = if_else(gene_symbol == "P", "P/V/C", gene_symbol))%>%
mutate(Viral_epitope_tile = paste0 (gene_symbol, pos_start, "_", pos_end )) %>%
filter (hfc > hfc_value)
Epitopes_list <- unique(Epitopes_list$Viral_epitope_tile)
Epitopes_dataset<- datatable %>%
filter(taxon_species == virus) %>%
mutate(Viral_epitope_tile = paste0 (gene_symbol, pos_start, "_", pos_end ))%>%
filter(Viral_epitope_tile %in% Epitopes_list) %>%
group_by(Viral_epitope_tile, DPI, Treatment)%>%
mutate(averaged_hfc = mean (hfc))
return (Epitopes_dataset)
}
Epitopes_per_group_graph_function <- function (datatable){
ggplot(datatable, aes(x = DPI , y = Viral_epitope_tile))+
geom_raster(aes(fill = averaged_hfc)) +
scale_fill_gradientn(colors= c("white", "darkblue", "black")) +
labs(x = "DPI") +
facet_wrap(~Treatment, ncol = 4, nrow = 1, scales = "fixed") +
theme_minimal() +
theme(axis.text = element_text(size=10))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
}
Epitopes_per_RM_per_individual_group_graph_function <- function (datatable, treatment_type ){
ggplot(subset(datatable, Treatment %in% treatment_type), aes(x = DPI , y = Viral_epitope_tile))+
geom_raster(aes(fill = hfc)) +
scale_fill_gradientn(colors= c("white", "darkblue", "black")) +
labs(x = "DPI") +
facet_wrap(~Treatment, ncol = 4, nrow = 1, scales = "fixed") +
theme_minimal() +
theme(axis.text = element_text(size=10))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
}
Epitopes_per_RM_per_group_graph_function <- function (datatable, treatment_type ){
ggplot(subset(datatable, Treatment %in% treatment_type), aes(x = DPI , y = Viral_epitope_tile))+
geom_raster(aes(fill = hfc)) +
scale_fill_gradientn(colors= c("white", "darkblue", "black")) +
labs(x = "DPI") +
facet_manual(Treatment~Macaques_ID, design=design) +
theme_minimal() +
theme(axis.text = element_text(size=10))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
}
Measles_epitopes_dataset <- Epitopes_list_and_dataset_function (Merged_hfc_all, "Measles virus", 2)
Epitopes_per_group_graph_function(Measles_epitopes_dataset)
Epitopes_per_RM_per_individual_group_graph_function(Measles_epitopes_dataset,"WTMeV")
Epitopes_per_RM_per_individual_group_graph_function(Measles_epitopes_dataset,"RDV_PEP")
Epitopes_per_RM_per_individual_group_graph_function(Measles_epitopes_dataset,"RDV_Late")
Epitopes_per_RM_per_individual_group_graph_function(Measles_epitopes_dataset,"Sham_control")
#Measles_epitopes_list_hfc_over_50
Measles_epitopes_dataset_hfc_over_50 <- Epitopes_list_and_dataset_function (Merged_hfc_all, "Measles virus", 50)
Epitopes_per_group_graph_function(Measles_epitopes_dataset_hfc_over_50)
Epitopes_per_RM_per_group_graph_function(Measles_epitopes_dataset_hfc_over_50,RDV_groups)
#Measles_epitopes_dataset_filter_H
Measles_epitopes_dataset_filter_H <- Measles_epitopes_dataset %>%
filter(Viral_epitope_tile != "P/V449_504" & Viral_epitope_tile != "N477_525" & Viral_epitope_tile != "H1_56") %>%
filter(grepl('H', Viral_epitope_tile))
Epitopes_per_group_graph_function(Measles_epitopes_dataset_filter_H)
Epitopes_per_RM_per_group_graph_function(Measles_epitopes_dataset_filter_H,RDV_groups)
#Measles_epitopes_dataset_filter_F
Measles_epitopes_dataset_filter_F <- Measles_epitopes_dataset %>%
filter(Viral_epitope_tile != "P/V449_504" & Viral_epitope_tile != "N477_525" & Viral_epitope_tile != "H1_56") %>%
filter(grepl('F', Viral_epitope_tile))
Epitopes_per_group_graph_function(Measles_epitopes_dataset_filter_F)
Epitopes_per_RM_per_group_graph_function(Measles_epitopes_dataset_filter_F,RDV_groups)
Epitope_tile_hfc_function <- function (epitope_dataset, virus, gene_symbol_variable, viral_epitope_tile_variable) {
epitope_tile_hfc_variable<- epitope_dataset %>%
filter (gene_symbol ==gene_symbol_variable) %>%
filter (Viral_epitope_tile ==viral_epitope_tile_variable) %>%
filter (Duplicate != 2) %>%
group_by(DPI, Macaques_ID, Treatment) %>%
mutate(sum_of_hfc = sum(hfc))
plot_list <<- list()
plot1 <- ggplot(epitope_tile_hfc_variable, aes(x=DPI, y=sum_of_hfc, group = Macaques_ID, color = Treatment)) +
geom_point() +
geom_line(linetype = "solid", color = "gray") + theme_bw() +
labs(x= "DPI" , y="peptide_hfc" , title = paste0("Hit-fold change of ", virus," ", gene_symbol_variable," protein\n", "peptide tile " , viral_epitope_tile_variable))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
facet_wrap(~Treatment)
plot2 <- ggplot(epitope_tile_hfc_variable, aes(x=DPI, y=sum_of_hfc, color = Treatment)) +
geom_point(position=position_jitterdodge(jitter.width = 0, jitter.height = 0))+
geom_boxplot()+
theme_bw() +
labs(x= "DPI" , y="peptide_hfc" , title = paste0("Average hit-fold change of ", virus," ", gene_symbol_variable," protein\n", "peptide tile " , viral_epitope_tile_variable))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
facet_wrap(~Treatment)
plot3 <- ggplot(epitope_tile_hfc_variable, aes(x=DPI, y=sum_of_hfc, group = Macaques_ID, color = Treatment)) +
geom_point() +
geom_line(linetype = "solid", color = "gray") + theme_bw() +
labs(x= "DPI" , y="peptide_hfc" , title = paste0("Hit-fold change of ", virus," ", gene_symbol_variable," protein\n", "peptide tile " , viral_epitope_tile_variable))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
facet_manual(Treatment~Macaques_ID, design=design)
plot_list <<-list(plot1, plot2, plot3)
}
Epitope_tile_hfc_function (Measles_epitopes_dataset, "measles", "N", "N477_525")
grid.arrange(grobs = plot_list, widths = c(2, 1, 1),layout_matrix = rbind(c(1, 3, 3), c(2, 3, 3)))
Epitope_tile_hfc_function (Measles_epitopes_dataset, "measles", "H", "H1_56")
grid.arrange(grobs = plot_list, widths = c(2, 1, 1),layout_matrix = rbind(c(1, 3, 3), c(2, 3, 3)))
Epitope_tile_hfc_function (Measles_epitopes_dataset, "measles", "P/V", "P/V449_504")
grid.arrange(grobs = plot_list, widths = c(2, 1, 1),layout_matrix = rbind(c(1, 3, 3), c(2, 3, 3)))
# Top viral species that overlapped in both VarScore and PPH dataset
overlapping_viral_species <- intersect(top_viral_species_human_lib_2, top_viral_species_human_lib_10)
overlapping_viral_species
## [1] "BK polyomavirus" "Hepatitis delta virus"
## [3] "Hepatitis delta virus genotype I" "Human respiratory syncytial virus"
## [5] "JC polyomavirus" "Measles virus"
## [7] "Simian virus 40"
Average_peptide_hits_per_protein_of_selected_viral_species_function <- function (viral_species_selected) {
hfc_dataset_of_selected_viral_species <- Viral_peptide_hits_datatable_function(Merged_hfc_all, viral_species_selected)
Average_peptide_hits_per_protein_function(hfc_dataset_of_selected_viral_species, viral_species_selected, "springgreen4")
}
Average_peptide_hits_per_protein_of_selected_viral_species_function("BK polyomavirus")
Average_peptide_hits_per_protein_of_selected_viral_species_function("Hepatitis delta virus")
Average_peptide_hits_per_protein_of_selected_viral_species_function("Human respiratory syncytial virus")
Average_peptide_hits_per_protein_of_selected_viral_species_function("JC polyomavirus")
Average_peptide_hits_per_protein_of_selected_viral_species_function("Simian virus 40")
Average_peptide_hits_per_protein_of_selected_viral_species_function("Simian foamy virus type 1")
Average_peptide_hits_per_protein_of_selected_viral_species_function("Simian foamy virus type 3")